rm(list=ls())

library(tidyverse)
library(readr)
library(stargazer)
library(stringr)
library(readxl)
library(texreg)
library(broom)
library(fastDummies)
library(xtable)
library(grid)
library(mediation)
library(gtable)
library(conflicted)
library(estimatr)
library(broom)
library(sjPlot)
library(gridExtra)
library(ggplot2)
library(ggpubr)

conflict_prefer("filter", "dplyr")
conflict_prefer("select", "dplyr")

dat <- read_csv("tass-wave-1-october-2023.csv")
dat <- dat %>%
  mutate(party_binary = ifelse(pid7 %in% c(1,2,3), "Democrats", 
                               ifelse(pid7 %in% c(5,6,7), "Republicans",
                                      ifelse(pid7 == 4, "Independent", "Other"))),
         source = recode(vignette_randomization,`1` = "N", `2` = "D", `3` = "R"),
         source = factor(source, levels = c("N", "D", "R")),
         party_binary = factor(party_binary, levels = c("Democrats", "Republicans", "Independent", "Other")))

######
## H1: Ban social media platforms
## H2: inpect social media content
## H3: prosecute foreign software
## H4: Foriegn government ties
## H5: President and the supreme court
## H6: limit immigration from China

dat$Q1 <- dat$H1
dat$Q2 <- dat$H2
dat$Q3 <- dat$H3
dat$Q4 <- dat$H4
dat$Q5 <- dat$H5
dat$Q6 <- dat$H6

summ <- dat %>%
  group_by(source) %>%
  summarise(across(Q1:Q6, ~mean(., na.rm=T), .names = "mean_{.col}"),
            across(Q1:Q6, ~sd(., na.rm=T), .names = "sd_{.col}"))

summ_longer_mean <- summ %>%
  select(source,mean_Q1:mean_Q6) %>%
  pivot_longer(cols = (-source),
               names_to = "DV",
               names_prefix = "mean_",
               values_to = "mean")

summ_longer_sd <- summ %>%
  select(source,sd_Q1:sd_Q6) %>%
  pivot_longer(cols = (-source),
               names_to = "DV",
               names_prefix = "sd_",
               values_to = "sd")

summ_longer <- bind_cols(summ_longer_mean, sd=summ_longer_sd$sd)
summ_longer <- summ_longer %>%
  arrange(DV)

summ_longer %>%
  ggplot() +
  geom_bar(aes(x=source, y=mean), stat="identity") +
  facet_wrap(~DV)

tr_effects <- tidy(lm_robust(cbind(Q1,Q2,Q3,Q4,Q5,Q6) ~ source, data = dat))
tr_effects %>%
  filter(term %in% c("sourceD", "sourceR")) %>%
  ggplot(aes(factor(term,labels=c('D','R')),
             estimate))+
  geom_point(size = 2,position=position_dodge(width=0.5),color='#A51417') +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.05, linewidth = 1.2,
                position=position_dodge(width=0.5),color='#A51417')+
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_y_continuous(breaks = -0.5:0.5) +
  ylim(-0.5:0.5) +
  facet_wrap(~outcome)+
  theme_bw()+
  labs(x='Source',y='Effect')+
  theme(strip.text = element_text(size=20),
        axis.text.y = element_text(size=20),
        axis.text.x = element_text(size=20),
        axis.title=element_text(size=20),
        legend.text=element_text(size=20),
        legend.title = element_text(size=20),
        title = element_text(size=30),
        panel.border = element_rect(color = "black", fill = NA, size = 1))

ggsave("Plots/001_main_effects.pdf", height = 8, width = 12)

#### by party
summ_byparty <- dat %>%
  group_by(source, party_binary) %>%
  summarise(across(Q1:Q6, ~mean(., na.rm=T), .names = "mean_{.col}"),
            across(Q1:Q6, ~sd(., na.rm=T), .names = "sd_{.col}"),
            across(Q1:Q6, ~sum(!is.na(.)), .names = "n_{.col}"),)

summ_byparty_longer_mean <- summ_byparty %>%
  select(source,party_binary,mean_Q1:mean_Q6) %>%
  pivot_longer(cols = (-c(source,party_binary)),
               names_to = "DV",
               names_prefix = "mean_",
               values_to = "mean")

summ_byparty_longer_sd <- summ_byparty %>%
  select(source,party_binary,sd_Q1:sd_Q6) %>%  
  pivot_longer(cols = (-c(source,party_binary)),
               names_to = "DV",
               names_prefix = "sd_",
               values_to = "sd")

summ_byparty_longer_n <- summ_byparty %>%
  select(source,party_binary,n_Q1:n_Q6) %>%  
  pivot_longer(cols = (-c(source,party_binary)),
               names_to = "DV",
               names_prefix = "n_",
               values_to = "n")

summ_byparty_longer <- bind_cols(summ_byparty_longer_mean, 
                                 sd=summ_byparty_longer_sd$sd,
                                 n=summ_byparty_longer_n$n)

summ_byparty_longer <- summ_byparty_longer %>%
  arrange(DV)

summ_byparty_longer_DR <- summ_byparty_longer %>%
  filter(party_binary %in% c("Democrats", "Republicans"))

summ_byparty_longer_DR %>%
  ggplot() +
  geom_bar(aes(x=source, y=mean, fill=party_binary), 
           stat="identity", position=position_dodge()) +
  ylim(0,6) +
  scale_fill_manual(values=c('#00AEF3','#E81B23')) +
  facet_wrap(~DV)

tr_effects_byparty <- tidy(lm_robust(cbind(Q1,Q2,Q3,Q4,Q5,Q6) ~ source*party_binary, 
                                     data = dat %>% filter(party_binary %in% c("Democrats", "Republicans"))))

## H1: ban social media platforms
h1plot <- plot_model(lm_robust(Q1 ~ source*party_binary,
                     data = dat %>% filter(party_binary %in% c("Democrats", "Republicans"))), 
           type = "int", terms = c("sourceD", "party_binary"), title = "Q1: ban social \n media platforms", colors = c("#00AEF3", "#E81B23"),
           axis.title = c("Source Partisanship", "Q1 Response"), legend.title = "Respondent \n Partisanship",
           show.legend = F)
# ggsave("Plots/011_her_effects_h1.pdf", height = 4, width = 5)

## H2: inspect social media content
h2plot <- plot_model(lm_robust(Q2 ~ source*party_binary,
                     data = dat %>% filter(party_binary %in% c("Democrats", "Republicans"))), 
           type = "int", terms = c("sourceD", "party_binary"), title = "Q2: inspect social \n media content", colors = c("#00AEF3", "#E81B23"),
           axis.title = c("Source Partisanship", "Q2 Response"), legend.title = "Respondent \n Partisanship",
           show.legend = F)
# ggsave("Plots/011_her_effects_h2.pdf", height = 4, width = 5)


## H3: prosecute foreign software
h3plot <- plot_model(lm_robust(Q3 ~ source*party_binary,
                     data = dat %>% filter(party_binary %in% c("Democrats", "Republicans"))), 
           type = "int", terms = c("sourceD", "party_binary"), title = "Q3: prosecute foreign \n software usage", colors = c("#00AEF3", "#E81B23"),
           axis.title = c("Source Partisanship", "Q3 Response"), legend.title = "Respondent \n Partisanship",
           show.legend = F)
# ggsave("Plots/011_her_effects_h3.pdf", height = 4, width = 5)

## H4: Foriegn government ties
h4plot <- plot_model(lm_robust(Q4 ~ source*party_binary,
                     data = dat %>% filter(party_binary %in% c("Democrats", "Republicans"))), 
           type = "int", terms = c("sourceD", "party_binary"), title = "Q4: foreign ties \n ineligible for office", colors = c("#00AEF3", "#E81B23"),
           axis.title = c("Source Partisanship", "Q4 Response"), legend.title = "Respondent \n Partisanship",
           show.legend = F)
# ggsave("Plots/011_her_effects_h4.pdf", height = 4, width = 5)

## H5: President and the supreme court (incumbent is Biden)
h5plot <- plot_model(lm_robust(Q5 ~ source*party_binary,
                     data = dat %>% filter(party_binary %in% c("Democrats", "Republicans"))), 
           type = "int", terms = c("sourceD", "party_binary"), title = "Q5: Biden can \n override SC", colors = c("#00AEF3", "#E81B23"),
           axis.title = c("Source Partisanship", "Q5 Response"), legend.title = "Respondent \n Partisanship",
           show.legend = F)
# ggsave("Plots/011_her_effects_h5.pdf", height = 4, width = 5)

## H6: limit immigration from China
h6plot <- plot_model(lm_robust(Q6 ~ source*party_binary,
                     data = dat %>% filter(party_binary %in% c("Democrats", "Republicans"))), 
           type = "int", terms = c("sourceD", "party_binary"), title = "Q6: limit immigration \n from China", colors = c("#00AEF3", "#E81B23"),
           axis.title = c("Source Partisanship", "Q6 Response"), legend.title = "Respondent \n Partisanship",
           show.legend = F)
# ggsave("Plots/011_her_effects_h6.pdf", height = 4, width = 5)

ggarrange(h1plot, h2plot, h3plot, h4plot, h5plot, h6plot,
                    ncol = 2, nrow = 3)

ggsave("Plots/011_her_effects.pdf", height = 8.4, width = 6)
